public abstract class DifferentialEquationSolver implements IDifferentialEquationSolver {
    protected IOdeRhs yPrime;
    protected double[] y0;
    protected double t0;
    
    // TODO:  Assumes y0, result of yPrime same size.  Add a check, throw an exception if not.
    public DifferentialEquationSolver(IOdeRhs yPrime, double[] y0, double t0) {
        int i;
    	this.yPrime = yPrime;
    	this.y0 = y0.clone();
    	this.t0 = t0;
    }
    
    public static boolean isFinite ( double[] v ) {
        int i ;
        for ( i = 0 ; i < v.length ; i ++ ) {
            if ( Double.isInfinite ( v[i] ) || Double.isNaN( v[i] ) ) {
                return false;
            }
        }
        return true;
    }
    
    // protected because assumes v1 and v2 are the same length
    protected static double[] add(double[] v1, double[] v2) {
        double[] result = new double[v1.length];
        int i;
        for ( i = 0; i < v1.length; i++ ) {
            result[i] = v1[i] + v2[i] ;
        }
        return result;
    }

    // protected because assumes v1 and v2 are the same length
    protected static double[] subtract(double[] v1, double[] v2) {
        double[] result = new double[v1.length];
        int i;
        for ( i = 0; i < v1.length; i++ ) {
            result[i] = v1[i] - v2[i] ;
        }
        return result;
    }    
    
    // calculate the 2-norm of a vector v
    public static double norm(double[] v) {
        double sum = 0;
        int i;
        for (i=0; i<v.length; i++) {
            sum += v[i] * v[i];
        }
        return Math.sqrt(sum);
    }

    public double[] advanceTo(double t, double h) {
        double dt = t - t0;
        int n = (int) Math.round(dt/h);
        int i;
        h = dt/n;
        for (i=0; i<n; i++) {
            step(h);
        }
        return y0.clone();
    }
    
    protected static double[] multiply(double c, double[] v) {
        double[] result = new double[v.length];
        int i;
        for ( i = 0 ; i < v.length; i++ ) {
            result[i] = c * v[i];
        }
        return result;
    }
}
